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Abstract 

An  anisotropic  phase-field  model  is  used  to  calculate  numerically  the  solidification 
patterns  of  a pure  material  into  an  undercooled  liquid  in  a two-dimensional  rectangular 
region.  In  the  phase-held  approach,  the  solid-liquid  interface  is  treated  as  diffuse,  and 
a dynamic  equation  for  the  phase  variable  is  introduced  in  addition  to  the  equation  for 
heat  flow.  The  phase-held  model  equations  are  solved  using  finite-difference  techniques 
on  a uniform  mesh.  Calculations  for  dendritic  growth  are  presented  for  both  four-fold 
and  six-fold  anisotropy,  and  the  effect  of  the  level  of  anisotropy  on  the  growth  of  a 
dendrite  is  investigated.  A previous  study  has  shown  that  performing  computations 
with  an  interface  that  is  sufficiently  thin  for  the  numerical  solution  to  accurately  rep- 
resent a sharp  interface  model  is  computationally  demanding.  However,  even  with  a 
relatively  thick  interface,  the  computations  using  the  phase-held  model  show  many  of 
the  qualitative  features  of  dendritic  growth,  and  the  method  is  well  suited  for  handling 
the  evolution  of  very  complex,  realistic  interface  shapes. 
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1.  Introduction 


Phase-field  models  provide  a comparatively  new  setting  for  the  theoretical  investigation  of 
phase  transitions,  such  as  solidification  of  a pure  material.  The  classical  free  boundary  ap- 
proach models  phase  change  between  two  phases  separated  by  a moving  interface  upon  which 
boundary  conditions  must  be  prescribed.  In  contrast,  the  phase- field  method  introduces  a 
continuous  transition  between  the  two  phases  across  a thin  layer  of  finite  thickness.  An  addi- 
tional variable  called  the  phase  field  is  introduced  whose  value  at  a point  identifies  the  phase. 
The  advantage  of  this  approach  lies  in  its  coherent  description  of  phase  change,  whereby  no 
distinction  is  made  in  the  model  between  the  solid,  liquid  or  interfacial  layer.  An  important 
consequence  of  the  phase-field  approach  is  that  numerical  solutions  of  the  resulting  equations 
can  yield  complicated  growth  morphologies  whose  connectedness  changes  in  time;  thus,  it  is 
well  suited  for  the  numerical  simulation  of  complicated  phenomena  such  as  dendritic  growth, 
a situation  which  is  very  difficult  to  solve  numerically  employing  algorithms  based  on  the 
classical  free-boundary  problem  formulation. 

Early  work  on  the  development  of  phase-field  models  for  first-order  phase  transitions 
can  be  traced  back  to  a number  of  investigators  [1-5].  The  development  of  a field  equation 
for  an  order  parameter  or  phase-field  variable  that  can  be  used  to  model  solidification  bor- 
rows ideas  from  the  related  physical  problems  of  diffuse  phase  boundaries  and  critical  point 
phenomena  [6-9].  Several  computations  using  a phase-field  model  have  been  performed 
for  one- dimensional  planar  or  radially  symmetric  geometries  [10-13].  Recently,  Kobayashi 
[14, 15]  reported  computations  of  a phase-field  model  for  a pure  material  in  two  and  three 
spatial  dimensions,  which  clearly  showed  the  evolution  of  solid  dendritic  structures  into  an 
undercooled  melt.  Kobayashi ’s  work  was  the  first  qualitative  demonstration  of  the  possible 
utility  of  phase-field  models  of  solidification  as  a computational  tool  for  modeling  compli- 
cated, realistic  solid/liquid  interfaces. 

More  recently,  a phase-field  model  has  been  developed  and  used  to  numerically  simu- 
late the  solidification  of  a pure  material  [16,17].  This  model  was  formulated  in  a manner 
consistent  with  irreversible  thermodynamics  [18]  and  careful  attention  was  paid  to  the  deter- 
mination of  the  phase-field  model  parameters  in  terms  of  the  classical  material  parameters. 
This  recent  model  is  also  well  suited  to  numerical  computation.  There  has  also  been  some 
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recent  work  on  the  development  of  a phase-field  model  for  the  solidification  of  an  alloy  in  an 
isothermal  system  [19,20].  The  methodology  employed  in  developing  these  newer  models  has 
provided  the  necessary  framework  for  the  development  of  a complete  phase- field  description 
for  the  more  general  treatment  of  a nonisothermal  binary  alloy. 

In  this  paper,  additional  computations  for  a two-dimensional  region  are  performed  using 
the  model  developed  recently  for  a pure  material.  The  time- dependent  heat  and  phase-field 
equations  are  solved  numerically  using  finite-difference  techniques  on  a uniform  computa- 
tional mesh  for  conditions  simulating  dendritic  growth  into  an  undercooled  liquid.  Compu- 
tations are  performed  for  both  four-fold  and  six-fold  anisotropy.  Some  qualitative  aspects  of 
using  the  phase-field  approach  for  such  calculations  is  investigated.  This  work  is  part  of  an 
ongoing  effort  to  develop  and  evaluate  phase-field  models  of  solidification. 


2.  Governing  Equations 


Phase-field  models  of  solidification  are  typically  derived  from  a Landau-Ginzburg  or  Cahn- 
Hilliard  type  of  free-energy  functional: 


T 


-L 


MT)  + -(vtf 


dV, 


(1) 


where  V is  the  region  occupied  by  the  system,  0(x,  t ) is  the  phase-field,  T(x,  t ) is  the  temper- 
ature and  e is  a gradient  energy  coefficient,  which  is  constant  for  an  isotropic  material.  The 
free-energy  density  /(<^>,  T)  is  a monotonic  function  of  T and  is  a double-well  with  respect 
to  c/. ),  which  yields  two  minima  with  respect  to  (f)  representing  the  free  energy  of  the  solid 
and  liquid  phases  for  each  temperature.  Different  choices  for  the  precise  functional  form  of 
/ have  been  used  in  the  various  earlier  models  [4,15].  In  most  of  the  previous  developments 
of  phase-field  models,  the  free  energy  Eq.  (1)  has  been  used  to  obtain  a kinetic  equation  for 
the  phase  field  by  requiring  that  it  evolves  in  a manner  such  that  the  total  free  energy  T 
decreases  monotonically  in  time.  The  simplest  form  of  a kinetic  equation  that  is  consistent 
with  this  requirement  is  the  following: 


d(j) 

dt 


6F 


(2) 


where  M\  is  a parameter  that  can  be  related  to  the  interface  kinetic  coefficient  of  the  sharp 
interface  problem.  Typically  in  the  formulation  of  a phase-field  model,  a modified  form  of 
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the  heat  equation  is  assumed  which  includes  an  appropriate  source  term  to  account  for  the 
liberation  of  latent  heat,  e.g., 


£Z!  -l  K—  — kV*t 
dt  +K  at  j’ 


(3) 


where  if  is  a constant  proportional  to  the  latent  heat  per  unit  volume,  and  k is  the  thermal 
diffusivity. 

Penrose  and  Fife  [18]  have  addressed  the  problem  of  deriving  the  set  of  equations  in 
a more  rigorous  fashion  by  employing  an  an  entropy  functional,  which  is  the  appropriate 
thermodynamic  potential  for  nonisothermal  situations.  They  derived  the  set  of  phase-field 
equations  using  a form  for  the  potential,  /(</>,  T),  introduced  in  earlier  studies  [3,4].  However, 
their  choice  for  the  form  of  the  potential  function,  / has  the  disadvantage  that  the  values 
of  <j)  representing  the  solid  and  liquid  states  depends  upon  temperature  and  other  model 
parameters.  This  limitation  complicates  the  implementation  of  the  phase-field  model  for 
practical  computation. 

A recently  developed  phase-field  model  is  used  for  the  calculations  here  which  is  derived 
in  [17].  The  model  is  derived  from  an  entropy  functional  in  a manner  similar  to  Penrose 
and  Fife  [18],  but  it  is  constructed  so  that  the  states  (f)  = 0 and  <f>  = 1 correspond  to  the 
bulk  solid  and  liquid  phases,  respectively,  independent  of  the  temperature.  This  has  the 
advantage  that  latent  heat  is  liberated  only  in  the  interfacial  region.  It  is  noted  here  that 
different  conventions  are  used  for  the  values  of  (f)  corresponding  to  the  bulk  phases  in  the 
various  phase-field  models.  This  new  phase- field  model  results  in  the  following  dimensionless 
governing  equations  for  the  temperature  and  phase  field: 


du  30  g{<f>)d<l>  2 

m + =Vu 

~mWt  = 3°s(4>)iaSu  - jg'(<t>)  + e2V2^.  (5) 

where  g((f>)  = <^>2(1  — (j))2  and  prime  denotes  differentiation.  Here,  length  has  been  scaled  on 
some  reference  length  scale  w of  say  the  dimensions  of  the  domain,  time  on  the  corresponding 
thermal  diffusion  time  w2  / n , and  temperature  by  putting  T — Tm  + A Tu,  where  Tm  is 
the  equilibrium  melting  temperature  and  AT  is  a reference  temperature  difference  (such 
as,  between  the  melting  temperature  and  the  initial  temperature  at  the  boundary  of  the 
domain). 
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This  isotropic  phase-field  model  is  characterized  by  four  dimensionless  parameters: 


c cAT 

L ’ 

(6) 

sflwL2 

(7) 
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K,  L 

8 

(9) 

e = — , 

where  c is  the  specific  heat  per  unit  volume,  L is  the  latent  heat  per  unit  volume,  a is  the  sur- 
face energy,  /i  is  the  interface  kinetic  coefficient,  and  8 is  a measure  of  the  interface  thickness, 
which  has  been  related  to  the  gradient  energy  coefficient,  e,  by  analyzing  a one-dimensional 
equilibrium  solution  of  the  dimensional  phase-field  equation.  In  [17],  it  is  shown  how  the 
phase-field  parameters  a and  m are  related  to  the  physical  parameters  which  characterize  the 
interface  dynamics  (i.e. , interfacial  energy  and  kinetic  coefficient).  Once  the  characteristic 
length  scale  w has  been  chosen,  knowledge  of  the  physical  properties  leaves  one  degree  of 
freedom,  namely  e,  which  then  is  used  to  set  the  interface  thickness.  It  is  expected  that  in 
order  to  model  the  physical  behavior  correctly,  the  interface  thickness  must  be  sufficiently 
small  compared  to  the  interfacial  macrostructures  that  are  to  be  modeled;  however,  from  a 
computational  viewpoint,  it  is  desirable  for  the  interface  thickness  to  be  as  large  as  possible 
in  order  that  accurate  solutions  of  the  phase-field  equations  can  be  obtained  for  practical 
computational  effort. 

In  order  to  simulate  dendritic  growth,  the  phase-field  equation  given  by  Eq.  (5)  is  modified 
to  include  anisotropy  in  the  parameter  e by  introducing 


e(6)  = er](9)  = e(l  + 7 cos  k6). 


The  angle  6 is  defined  as  the  angle  between  the  normal  to  the  interface  and  the  x-axis,  and  k 
specifies  the  mode  number.  Rederiving  Eq.  (5)  after  including  the  variation  of  the  gradient 
energy  coefficient  with  orientation  yields  a modified  version  of  the  phase-field  equation: 


(f>-^  + 30eaS(f){l-(t))u  - e2-^-  (v(8)v'{0)^ 


(*> W(^)  + ' fo2Wv«}  • (10) 
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In  the  phase-field  model,  the  interface  is  represented  by  a range  of  contour  levels  of  <f),  e.g., 
0-1  < (f)  < 0.9.  In  order  to  compute  the  anisotropic  behavior,  the  orientation  angle  is 
determined  in  terms  of  the  phase  field  4>  using  the  following  relation  for  the  normal  vector: 


n = 


V<£ 

m 


= cos  6x  + sin  9y. 


From  this  expression  we  have  the  definition 


tan  6 = 

f Vx 


and,  in  addition,  we  obtain 


Q _ fyxfyxy  fiyfixx  Q ( fix&yy  tyyfyxy 

x ~ m2  ’ v " m2  ■ 

These  relations  are  used  to  compute  the  terms  in  Eq.  (10)  which  arise  from  expanding  the 
derivatives  of  rj(6). 

By  conducting  an  asymptotic  analysis  of  the  phase-field  equations  in  the  limit  e — * 0 with 
5,  a,  and  m of  order  one,  the  classical  description  of  an  appropriate  free  boundary  problem 
is  recovered.  The  details  are  analogous  to  those  given  by  Caginalp  [5]  for  the  form  of  the 
free  energy  investigated  by  him,  and  a thorough  development  of  the  interfacial  conditions 
for  the  present  phase-field  model  with  the  assumed  form  for  e(9)  is  given  in  [21].  The  free 
boundary  problem  satisfied  by  the  leading-order  temperature  for  the  phase-field  model  given 
by  Eq.  (4)  and  Eq.  (10)  is 


(ii) 


with  interfacial  boundary  conditions 


du 

dn 


Liquid 

Solid 


(12) 


and 

“ = _r{^  + I,w  + ,w:}’  (13) 
where  T = aTm/[wLAT]  is  the  dimensionless  capillary  length,  vn  is  the  dimensionless  normal 

velocity  of  the  interface  into  the  liquid  and  /C  is  the  dimensionless  interfacial  curvature 

(measured  positive  for  convex  projections  into  the  liquid).  The  interfacial  boundary  condition 

Eq.  (13)  reduces  to  the  anisotropic  form  of  the  Gibbs-Thomson  equation  when  vn  — 0;  the 

term  including  vn  incorporates  the  effect  of  interface  kinetics.  Note  that  allowing  e to  depend 

on  6 modifies  both  the  term  proportional  to  velocity  and  the  term  proportional  to  curvature. 
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3.  Numerical  Computations 

The  objective  here  is  to  evaluate  the  behavior  of  the  phase-field  model  by  simulating  the 
growth  of  a two-dimensional  dendrite  into  an  undercooled  liquid.  In  order  to  perform  the 
simulation,  the  governing  equations  are  solved  numerically  in  a two-dimensional  rectangular 
domain  of  sufficient  size  to  allow  for  the  development  of  characteristic  dendritic  structure 
(e.g.,  side  arms).  Even  with  the  phase-field  approach,  an  optimum  solution  procedure  should 
employ  some  type  of  adaptive  technique,  particularly  for  the  phase  variable,  since  it  varies 
over  a very  small  part  of  the  domain.  However,  we  chose  at  present  to  use  straightforward 
finite  difference  solution  techniques  applied  on  a uniform  computational  grid,  since  the  nu- 
merical implementation  of  such  techniques  is  straightforward  and  they  are  well  suited  for 
taking  full  advantage  of  highly  vectorized  large-scale  computers. 

The  governing  equations  given  by  Eq.  (4)  and  Eq.  (10)  are  a pair  of  coupled,  second- 
order,  nonlinear  parabolic  equations.  They  are  discretized  spatially  using  second-order  finite 
differences  on  a uniform  grid  characterized  by  mesh  spacings  AX  and  AY  in  the  x and  y 
coordinate  directions,  respectively;  for  the  temporal  discretization  we  introduce  the  time  step 
At.  In  order  to  maximize  the  computational  efficiency,  we  employ  explicit  time-differencing 
on  the  (f)  equation  which  is  nonlinear  in  all  terms  except  for  the  highest-order  spatial  deriva- 
tives. The  heat  Eq.  (4)  is  linear  in  the  temperature  u but  contains  the  source  term  depending 
on  (f).  For  the  present  calculations,  the  thermal  field  diffuses  more  rapidly,  so  to  avoid  the 
diffusive  stability  requirement,  the  alternating-direction  implicit  method  (ADI)  is  used  on 
Eq.  (4).  The  methods  employed  here  are  described  in  many  standard  texts  on  finite  difference 
techniques  for  partial  differential  equations  (e.g.,  [22]). 

For  the  two-dimensional  rectangular  computational  domain,  vanishing  Neumann  condi- 
tions for  both  (j)  and  u are  applied  at  the  boundaries,  which  have  the  scaled  dimensions  of 
Xl  and  Yl  in  the  x and  y coordinate  directions,  respectively.  We  note  that  some  of  the 
results  presented  here  for  the  temperature  and  phase  field  are  displayed  after  reflecting  the 
computational  domain  about  the  line  y — 0,  which  corresponds  to  the  axis  of  the  dendrite. 
From  the  chosen  definition  of  the  dimensionless  variables,  the  value  u = 0 corresponds  to 
the  melting  temperature  of  the  pure  material,  while  u = — 1 is  the  initial  undercooling  tem- 
perature (the  dimensional  undercooling  is  AT).  Initially  for  each  calculation,  a small  region 


-7- 


of  solid  ((f)  — 0,u  — 0)  is  located  at  the  x — 0,  y = 0 corner  of  the  domain,  and  the  remain- 
der of  the  domain  is  undercooled  liquid  (u  = —1).  The  shape  of  the  initial  solid  region  is 
one-quarter  of  a circle  with  an  initial  radius  denoted  by  ra. 

The  two-dimensional  simulations  were  performed  using  property  values  for  pure  nickel; 
the  values  for  the  required  material  parameters  are  readily  available  (see  for  example  [20]). 
The  property  values  used  to  determine  the  dimensionless  parameters  for  the  computations  are 
a — 3.7xl0-5  J/ cm2,  Tm  = 1728  K,  L = 2350  J/ cm3,  c = 5.42  J/  K cm3,  k - 0.155  cm2/s, 
and  fi  — 285  cm/K  s.  In  order  to  completely  determine  the  dimensionless  parameters  Eq.  (7) 
- Eq.  (9),  we  must  choose  values  for  the  reference  length,  w,  and  the  interface  thickness  8. 
The  choice  of  these  parameters  is  based  on  the  physical  structure  we  wish  to  compute  and 
the  practical  limitations  of  accurately  resolving  gradients  within  the  interfacial  region  for  a 
desired  computational  domain  of  size  Xi  and  Yl- 

Because  the  computational  domain  is  thermally  insulated,  the  value  of  the  dimensionless 
undercooling  S determines  the  portion  of  the  domain  which  can  be  solidified.  The  computa- 
tions are  facilitated  by  larger  values  of  S (large  values  of  undercooling,  AT)  because  the  solid 
encompasses  a larger  portion  of  the  computational  domain.  In  addition,  a larger  value  of  S 
produces  a thinner  thermal  boundary  layer  in  the  liquid  so  the  effect  of  the  adiabatic  outer 
boundaries  is  less  pronounced.  A value  of  S = 0.5  is  used  for  all  the  calculations  presented. 
This  value  corresponds  to  an  actual  undercooling  of  217  K,  which  is  an  attainable  level  of 
undercooling  for  nickel.  Based  on  some  preliminary  computational  experiments  to  determine 
how  the  interface  thickness  and  the  size  of  the  domain  should  be  chosen  in  relation  to  the 
tip  radius,  we  chose  a length  scale  w — 2.1  x 10-4  cm  which  yields  the  parameter  value 
a = 400  in  the  definition  Eq.  (7).  For  the  physical  parameters  of  nickel  given  above,  the 
dimensionless  parameter  m has  the  value  0.05. 

With  the  parameters  S',  a and  m specified,  the  parameter  e determines  the  thickness  of 
the  interface.  The  computational  resolution  was  determined  by  the  extent  of  the  domain 
chosen  (Xl  and  Yl)  and  the  number  of  grid  points  used  in  the  discretization  of  the  domain. 
For  all  the  computations  presented  here,  e = 0.005  is  used;  this  value  was  found  to  give  an 
optimal  interface  width  from  the  point  of  view  of  computational  efficiency  for  the  present 
parameters  and  computation  scheme.  A more  complete  discussion  about  the  choice  of  the 
model  parameters  and  the  effect  of  varying  the  interface  thickness  is  given  in  [16]. 
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4.  Results  and  Discussion 


In  Wheeler  et  al.  [16],  a quantitative  evaluation  of  the  phase-field  model  for  one-  and  two- 
dimensional  computations  was  performed  by  investigating  the  effects  of  interface  width  and 
spatial  and  temporal  numerical  resolution.  For  one-dimensional  spherically  symmetric  calcu- 
lations, very  good  agreement  was  obtained  between  the  phase-field  results  and  computations 
from  the  corresponding  sharp  interface  formulation  solved  numerically  using  a Green’s  func- 
tion approach.  The  phase-field  model  was  evaluated  by  comparison  to  existing  sharp  interface 
theories  of  two-dimensional  dendritic  growth.  Exact  comparisons  were  not  possible  owing  to 
the  differences  in  the  details  included  in  the  theories,  but  the  phase-field  calculations  were 
shown  to  produce  behavior  which  is  in  agreement  with  the  more  well-established  theories, 
i.e.,  the  Ivantsov/microscopic  solvability  theories. 

Here,  we  present  some  additional  computations  which  display  the  complexity  of  the  so- 
lidification structures  that  can  be  produced  by  the  phase-field  model  and  show  that  they 
give  a qualitative  representation  of  dendritic  growth  behavior.  Similar  computations  were 
performed  by  Kobayashi  [15]  using  a related  phase-field  model.  The  contribution  of  the 
present  calculations  is  that  they  are  obtained  from  a phase-field  model  which  has  been  ana- 
lyzed in  considerable  detail  and  for  which  the  model  parameters  have  been  well-established 
[16,17,21], 

The  results  presented  are  displayed  as  contours  of  the  temperature  and  phase  field.  In 
all  cases,  for  each  temperature  plot  eleven  isotherms  are  shown  which  are  fractions  of  the 
maximum  magnitude  of  the  temperature  (0.05,  0. 1-0.9,  and  0.95).  In  the  phase  field  plots,  the 
interface  is  represented  by  contours  in  the  range  0.1  < (f)  < 0.9  which  are  indistinguishable. 

In  figs.  1 and  2,  we  present  computations  which  illustrate  the  role  played  by  anisotropy 
in  the  parameter  e.  Both  calculations  were  performed  on  half  the  domain  shown  with  Xl  = 
4.5  and  Yl  = 2.25  with  600  grid  points  in  the  x direction  and  300  points  in  the  y direction 
(AX  = AY  = 0.0075).  For  these  calculations,  a time-step  of  5 x 10-5  was  used.  In  the 
figures  for  the  purpose  of  illustration,  the  x direction  has  been  oriented  so  that  it  is  vertical 
and  the  solution  has  been  reflected  about  the  center  vertical  line  ( y = 0).  Fig.  1 shows  the 
temperature  and  phase-field  contours  at  three  dimensionless  time  levels  in  the  computation 
(t  = 0.1,  0.3,  and  0.5)  for  the  case  when  there  is  no  anisotropy  (7  = 0).  The  initially  circular 
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solid  region  is  subject  to  morphological  instability  and  breaks  down  into  a number  of  finger- 
like structures  which  grow  and  experience  further  instability  in  the  form  of  tip-splitting. 
With  four-fold  anisotropy  ( k = 4)  and  an  anisotropy  parameter  value,  7 = 0.02,  fig.  2 shows 
the  temperature  and  phase-field  contours  at  the  same  three  dimensionless  time  levels  as  in 
fig.  1.  From  the  definition  of  the  orientation  angle,  the  preferred  growth  directions  for  k = 4 
are  along  the  coordinate  axes.  Additionally,  there  is  a short  solid  finger  with  an  orientation 
of  approximately  45  degrees  which  has  grown  out  early  in  the  calculation.  After  an  initial 
transient  period,  the  tips  of  the  needle  crystals  reach  a well-defined  operating  state  with  a 
particular  radius  and  growth  velocity.  These  values  remain  constant  until  the  proximity  of 
the  outer  boundary  of  the  domain  inhibits  growth. 

In  fig.  3 isotherms  and  phase-field  contours  at  a single  time  level  of  0.6  are  plotted  for 
a lower  value  of  the  anisotropy  parameter,  7 = 0.01.  Again,  because  of  the  influence  of 
anisotropy,  a well-defined  needle  crystal  is  obtained.  However,  in  this  case,  regular  side- 
branching  occurs  at  a distance  behind  the  tip  of  the  growing  crystal.  The  occurrence  of  side 
branching  is  found  to  depend  on  the  anisotropy  level  (7  value),  and  to  some  extent  on  the 
numerical  error;  the  side-branched  structure  occurs  more  readily  as  the  discretization  error 
increases.  For  this  set  of  calculations,  the  mesh  spacing  AX  = AY  = 0.015  was  used  and  the 
time-step  was  1 x 10-4.  The  dependence  of  side-arm  formation  on  anisotropy  level  was  also 
displayed  by  the  simulations  of  Saito,  et  al.  [23],  which  were  performed  using  a boundary 
integral  approach  to  solve  a sharp  interface  model  with  a quasistationary  heat  equation  and 
no  kinetic  effects. 

Kobayashi  [15]  found  that  the  occurrence  of  side-branched  structure  is  enhanced  by  in- 
troducing random  noise  into  the  calculations.  For  the  subsequent  calculations  presented, 
random  noise  of  small  amplitude  has  been  introduced  into  the  temperature  term  in  the 
phase-field  equation  Eq.  (10),  which  represents  a disturbance  to  the  interface.  Fig.  4 is  a 
calculation  with  anisotropy,  7 = 0.01,  and  with  a 1%  random  noise  level.  Isotherms  and 
phase-field  contours  are  shown  at  four  time  levels  which  are  equally  spaced  except  for  the 
lowest  two  time  levels.  At  the  lowest  time  level  shown,  the  preferred  directions  along  the 
coordinate  axes  experience  the  fastest  growth  and  the  solid  no  longer  resembles  the  initially 
circular  region.  At  the  second  time  level,  side-arm  formation  has  begun  to  occur,  and  in  the 
last  two  time  levels  the  side-branched  dendritic  structure  is  clearly  evident.  The  presence 
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of  the  random  noise  does  not  effect  the  overall  growth  conditions  of  the  needle  crystal  (tip 
radius  and  velocity),  but  does  enhance  the  irregular  production  of  the  secondary  arms.  Be- 
cause the  actual  computational  domain  is  shorter  in  the  y direction,  the  growing  dendrite 
tip  reaches  the  outer  boundary  in  that  direction.  As  in  the  previous  calculations,  this  one 
was  performed  on  half  the  domain  shown  and  was  reflected  about  the  vertical  centerline. 
For  the  calculation,  Xl  — 9.0  and  Yl  = 4.5  with  600  grid  points  in  the  x direction  and  300 
points  in  the  y direction  yielding  AX  = AY  = 0.015.  The  time-step  was  1 x 10~4. 

Finally,  we  present  a set  of  calculations  with  six-fold  anisotropy  ( k = 6)  in  fig.  5.  Unlike 
the  previous  calculations,  a square  domain  with  Xl  = Yl  = 9 has  been  used.  There  are 
900  grid  points  in  each  direction  and  the  time-step  was  again  1 x 10-4.  The  calculation 
was  initiated  by  a small  circular  region  of  solid  located  in  the  center  of  the  domain.  The 
anisotropy  level  is  7 = 0.01,  and  a 1%  level  of  random  noise  is  used  to  facilitate  the  occurrence 
of  irregular  side  branching.  The  development  of  the  dendritic  structure  with  six-fold  primary 
structure  is  shown  in  the  figure  at  the  four  dimensionless  time  values:  0.1,  0.2,  0.4,  and  0.6. 
The  occurrence  of  secondary  side-branches  and  even  some  tertiary  branch  formation  is  clearly 
displayed  in  the  simulation  as  well  as  the  inclusion  of  a liquid  pocket.  Also  evident  in  this 
figure  and  in  the  previous  one  is  the  occurrence  of  dendrite  side-arm  coarsening.  All  these 
features  and  the  fact  that  the  calculations  produce  a well-defined  operating  state  for  the 
dendrite  tip  are  good  examples  of  how  the  phase-field  model  can  produce  realistic  dendritic 
growth  behavior. 

The  results  presented  here  demonstrate  the  complexity  of  the  solidification  morphology 
that  can  be  simulated  using  the  phase-field  approach.  This  work  is  part  of  an  ongoing  effort 
to  evaluate  the  behavior  of  phase-field  models  both  qualitatively  and  quantitatively,  in  order 
to  access  their  utility  as  a practical  tool  for  computing  solidification  growth  morphologies  of 
real  materials  under  a variety  of  conditions. 
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Figure  Captions 


Figure  1.  The  computed  isotherms  (left)  and  phase-field  contours  (right)  at  three  time 
levels:  0.1,  0.3,  and  0.5  (shown  from  bottom  to  top).  The  calculation  is  for  no  anisotropy,  7 
= 0.  Note  that  the  computations  are  symmetric  about  the  vertical  centerline. 

Figure  2.  The  computed  isotherms  (left)  and  phase-field  contours  (right)  at  three  time 
levels:  0.1,  0.3,  and  0.5  (shown  from  bottom  to  top).  The  calculation  has  four-fold  anisotropy 
( k = 4)  with  7 = 0.02.  Note  that  the  computations  are  symmetric  about  the  vertical 
centerline. 

Figure  3.  The  computed  isotherms  (left)  and  phase-field  contours  (right)  at  a single  time 
level  0.6  for  a calculation  with  four-fold  anisotropy  ( k = 4)  and  7 = 0.02.  Note  that  the 
computations  are  symmetric  about  the  vertical  centerline. 

Figure  4.  The  computed  isotherms  (left)  and  phase-field  contours  (right)  at  four  time  levels: 
0.1,  0.3,  0.9  and  1.5  (shown  from  bottom  to  top).  The  calculation  has  four-fold  anisotropy 
(k  = 4)  with  7 = 0.01  and  a 1%  random  noise  level.  The  computations  are  symmetric  about 
the  vertical  centerline. 

Figure  5.  The  computed  isotherms  (left)  and  phase-field  contours  (right)  at  four  time  levels: 
0.1,  0.2,  0.4  and  0.6  (shown  from  bottom  to  top).  The  calculation  has  six-fold  anisotropy  ( k 
— 6)  with  7 = 0.01  and  a 1%  random  noise  level. 
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